\vsssub
\subsubsection{~$S_{is}$: Floe-size dependent scattering and dissipation} \label{sec:IS2}
\vsssub

\opthead{IS2}{\ws}{F. Ardhuin, C. Sevigny, G. Boutin,  D. Dumont,  T. Williams}


The implementation of this scattering term generally 
follows the approach of \cite{art:MM06}, to which has been added an estimation of the breakup of the ice 
by waves to be able to update a maximum floe size diameter. Finally a creep-based dissipation was also 
combined with the scattering.  


The scattering source term is defined by the scattering coefficient $\beta_{\mathrm{is},\mathrm{MIZ}}(\theta-\theta')$ which by default is isotropic, but can be given any directional dependency by setting {\code IS2ISOSCAT=FALSE} in namelist {\F SIS2}. The source term is thus,
\begin{equation}
 \frac{S_{\mathrm{is}}(k,\theta)}{\sigma} =  \int_0^{2\pi}\beta_{\mathrm{is},\mathrm{MIZ}}(\theta-\theta') [s_\mathrm{scat}  N(k,\theta')-N(k,\theta)] d\theta'  .
 \end{equation}
where $s_\mathrm{scat}$ is set to 1.0 by default but can be modified by 
{\code IS2BACKSCAT} in namelist {\F SIS2}. \\


The determination of scattering coefficients  $\beta_{\mathrm{is},\mathrm{MIZ}}$ is based on the theoretical 
reflection coefficient $\alpha_n(\sigma,h)$
for waves with a normal incidence going from a half-plane of open water to a half-plane of ice-covered water with a 
constant ice thickness $h$.  
Values of $\alpha_n(\sigma,h)$ as computed by \cite{art:KM08}
 (if {\code IS2WIM1=0.}) or by \cite{art:Ben12} (if {\code IS2WIM1=1.})  are tabulated 
in the {\code W3SIS2MD} module. Following \cite{art:Dea11}, the broken ice is treated as a series of 
such ice-water interfaces. Neglecting multiple reflections, the scattering parameterization defines the attenuation 
per unit time as if the ice-covered part of a grid cell was a succession of floes of mean diameter $D_{m}$ with a partial 
reflection $\alpha_n(\sigma,h)$ for each floe, giving, 
\begin{equation}
\beta_{\mathrm{is},\mathrm{MIZ}}=c_i~c_g \alpha_n(\sigma,h) / D_{m},
 \end{equation}
where $c_i$ is the ice concentration. \\
 

The estimation of the mean floe diameter $D_m$ is based on an assumed power law for the number of 
floes of diameter $D$, taken  proportional to $D^{-\gamma}$. This power law is further assumed to 
apply for $D$ ranging from 
the minimum $D_{\min}$ and a maximum $D_{\max}$. The average is thus given by  
\begin{equation}
D_m=\frac{\gamma}{\gamma -1}\times\frac{D_{\max}^{-\gamma +1}-D_{\min}^{-\gamma +1}}{D_{\max}^{-\gamma}-D_{\min}^{-\gamma}}.
\label{analytic_Dbar}
\end{equation}

\noindent
At present $D_{\min}$ is a user-supplied value. $D_{\max}$ can either be 
provided as a forcing field, e.g. from an 
ice model or some observations, or, if the namelist parameter {\code IS2BREAK} is set to {\code TRUE}, 
estimated from 
the breaking of the ice by the local wave field. If the namelist parameter {\code IS2DUPATE}  is set to  
{\code TRUE}, small values of $D_{\max}$ will persist even if the waves become too small to be able to break the ice to 
that size. This is probably the proper model 
use when  external forcing/coupling is available (e.g. advection of ice properties in an ice model). 
 On the contrary, if {\code IS2DUPDATE} is set to  {\code FALSE}, the value of $D_{\max}$ will be always adjusted to the local sea state, 
even if that means increasing $D_{\max}$.\\


Ice breaking by waves of wavelength $\lambda$ is assumed to produce floes of diameter $\lambda / 2$. 
In the parametrization, ice breaking occurs if the three following criteria are fulfilled \citep{art:Wea13}:
\begin{enumerate}
\item $\lambda/2 \geq D_{\min} \quad \mathrm{and} \quad  \lambda/2 \leq D_{\max}$
\item $D_{\max}>D_c$, as it exists a critical diameter, which depends on ice properties, below which no flexural failure is possible
\item $\varepsilon>\varepsilon_c$, the strain due to the incoming wave has to 
be greater than a defined critical strain
\end{enumerate}
The first criterion is simply checking that the new value of $D_{\max}$ will be larger than $D_{\min}$ and smaller 
than the previous value of $D_{\max}$.\\

The second criterion relies on \cite{inc:M86} with the correction specified in \cite{art:Bea18}, where $D_{c}$ is defined as
\begin{equation}
D_{c}=\dfrac{1}{2}\left(\frac{\pi^4 Y^* h^3}{48 \rho g (1 -\nu ^2)}\right)^{1/4}. 
\end{equation}

The third criterion corresponds to the flexural strain threshold. The horizontal 
strain caused by waves is related to 
the curvature of the ice layer, which, in one dimension is 
$\varepsilon = 0.5 h \partial ^2 \eta_{\mathrm ice} / \partial x^2$. 
The strain variance is given 
\begin{equation}
\langle\varepsilon^2\rangle = \left( \frac{h}{2} \right) ^2   \int_{k_1}^{k_2} k_{ice}^4 F(k)dk ,
\label{strain_sign}
\end{equation}
where $h$ is the ice thickness and $k_{ice}$ is the wavenumber $2 \pi /\lambda_{ice}$. 
Borrowing from wave breaking 
ideas \citep{art:BBY00}, the integration of the curvature variance is limited around 
the local wavenumber $k_{\mathrm ice}$. 
We also note that we have defined an effective minimum ice thickness $h_{\min}$ so that, 
if $h<h_{\min}$, the strain 
variance is computed with $h=h_{\min}$ to avoid unbreakable elastic thin ice in the model 
that does not correspond to 
usual observations. We thus take $D_{\max}$ to be half the wavelength of the shortest 
waves for which the following criterion is met 
\begin{equation}
F_{break}\, \sqrt{\varepsilon^2} > \frac{\sigma_c}{Y^*} ,
\label{strain_crit}
\end{equation}
where $\sigma_c$ is the ice flexural strength. $F_{break}$ is a factor representing random waves and adjustable with 
the {\F SIS2} namelist parameter  {\code IS2BREAKF}. It should in theory depend on the duration for which the ice is forced by the 
waves, and, based on the typical maximum value over 500 Rayleigh-distributes waves, was taken to be $F_{break}=3.6$. 
$F_{break} \, E_s$ is thus the maximum strain for random waves. \\ 




Inelastic dissipation was added in this routine, following \cite{art:Wad73}, because it critically 
depends on the floe size.  It assumes that the floes deformation is not fully elastic, and that the secondary creep under the wave-induced cyclic causes the dissipation of wave energy into heat. 
We use the ice floe law
\begin{equation}
%\left(\frac{d\varepsilon}{dt}\right)_{ij}=\frac{\tau^{n-1}}{B^n}\sigma_{i,j}'
\left(\frac{d\varepsilon}{dt}\right)_{ij}=\frac{\tau^2}{B^3}\sigma_{i,j}',
\end{equation}
$B$ is the floe law  constant and is a function of ice temperature. Using the normalized parameter estimated by \cite{art:Cea98} from 
laboratory experiments, $A=10^{11}$,  and  a uniform ice temperature of 270~K  gives a value 
of $B=10^7$~s$^{1/3}$. The volumic dissipation rate is 
\begin{equation}
%\frac{de}{dt}= | \sigma_{xx}^{n+1}/(2B)^{n} |.
\frac{de}{dt}= | \sigma_{xx}^{4}/(2B)^3 |.
\end{equation}
Also, the cyclic deformation of the ice can require a much larger elastic energy than the gravity potential 
energy, but this is only true if the ice is not broken. As a result, working with a wave elevation spectrum $E(k)$ 
could introduce large changes in $E(k)$ when the ice is broken or reformed. Instead we prefer 
to work with an energy spectrum $ R C_g E(k)/C_{g,\mathrm{ice}}$, using the coefficient $R$ introduced by \cite{art:Wad73}, which 
is the ratio of elastic to gravity potential energies. 
%For monochromatic waves of amplitude $a$, the energy by unit surface of waves in ice is thus 
%$R= \frac{1}{2} \rho g a^2 C_g$. It results that, without reflection nor dissipation, the wave variance in ice is such as :
%\begin{equation}
%E_{ice}=\frac{E\,C_{g}}{Cg_{ice}\,R}
%\label{eq:RWad}
%\end{equation}
%Note that if creep is activated (hence if {\code IS2BREAKE} $ > 0$), eq.~\ref{eq:RWad} 
%is also applied for the scattering curvature variance computation in $S_{IS2}$. 
For unbroken ice  $R$ is 
\begin{equation}
R=1+C_R\frac{4Y^*h^3\pi^4}{3\rho g \lambda^4(1-\nu^2)},
\label{R}
\end{equation}
where we have been careful that \cite{art:Wad73} used $2h$ for the  ice thickness, 
and $C_R$ is by default set to 1.0 using the namelist parameter {\code IS2BREAKE}, but it can 
be set to zero to work with the true elevation spectrum instead. This factor $R$ is also 
applied in the calculation of ice breakup by the waves.\\


The inelastic dissipation is linearized as  $S_{\mathrm{ine}}=-\alpha_{\mathrm{ine}} E_{ice}$. 
The coefficient $\alpha_{\mathrm{creep}}$ was adapted from the 
\cite{art:Wad73} monochromatic formula and is equal to
\begin{equation}
%\alpha_{\mathrm{ine}}=0.05 B h^5 \left(\frac{Y^*}{2B(1-\nu^2)}\right)^{(n+1)} I_n k^{n+1} \frac{C_g^2}{\rho g C_{g_{ice}}R^2} F \int_{k_1}^{k_2} k_{ice}^4 E(k)dk 
\alpha_{\mathrm{ine}}=0.05 B h^5 \left(\frac{Y^*}{2B(1-\nu^2)}\right)^{4} I_3 k^4 \frac{C_g^2}{\rho g C_{g_{ice}}R^2} F_{\rm broken} \int_{k_1}^{k_2} k_{ice}^4 E(k)dk, 
\label{eq:alpha_creep}
\end{equation}
where $I_3=\frac{1}{\pi} \int_0^\pi \sin^{4}\beta d\beta$. Details of the computation are given in \cite{art:Bea18}.
$F_{\rm broken}$  is a heuristic smooth transition from unbroken to broken ice, so that the dissipation 
gradually goes to zero for waves much longer than the floe sizes, because in that case the ice does 
not deform and produces no dissipation of wave energy, 
\begin{equation}
F_{\rm broken}=\tanh \left( {\frac{D_{\max}-C_{\lambda}\lambda_{ice}} {D_{\max}C_{smooth} } }\right) .
\label{smooth}
\end{equation}
Inelastic dissipation is computed after updating $D_{\max}$. 
%dissipation is full for
%$\lambda_{ice}\leq D_{\max}$ and disappears for $C_{\lambda}\lambda_{ice}< D_{\max}$.  
The two parameters in this smooth transition $C_{\lambda}$ and $C_{smooth}$ are 
set to  $0.4$ and $0.2$ by the adjustable namelist parameters {\code IS2CREEPD} and {\code IS2CREEPC}.\\

Possibility of subsituting this inelastic dissipation by an anelastic dissipation term has also been added. It is activated by setting the namelist parameter {\code IS2ANDISB} to {\code TRUE}. 
Anelastic dissipation corresponds to the energy dissipated into heat during the oscillatory motion of the dislocations induced by the cyclic stress associated to waves.  This behaviour results in a hysteresis that can be seen in stress-strain diagrams when sea ice is submitted to a sinusoidal stress as presented in Fig.~4 of \cite{art:Cea98} study. This latter article also suggests a model for the anelastic behaviour of sea ice, with a stress-strain relationship that enables to compute the area within the ellipse which results from the hysteresis. Similarly to what has been done for inelastic dissipation, anelastic dissipation is linearized as  $S_{\mathrm{ane}}=-\alpha_{\mathrm{ane}} E_{ice}$. 
\cite{art:Bea18} have derived  $\alpha_{\mathrm{ane}}$  from \cite{art:Cea98} model for a monochromatic stress. They obtained the following coefficient: 
\begin{equation}
\alpha_{\mathrm{ane}}= \dfrac{A}{6}\left(k_i^{2}\dfrac{Y^*}{(1-\nu^{2})\rho gG}\right)^{2}h^{3}\frac{C_g}{GC_{g,i}} F_{\rm broken},
\label{eq:alpha_ane}
\end{equation}  
where  $F_{\rm broken}$ is the same as for inelastic dissipation and $A$ is equal to:
\begin{equation}
A=\dfrac{4}{3}\sigma \alpha_{d}~\delta D^{d}~\dfrac{1}{\exp(\alpha_{d}s)+\exp(-\alpha_{d}s)}, 
\label{eq:A_ane}
\end{equation}  
in which terms are detailed in the table at the end of this section.\\

Finally we recall the various model parameters used in IS2 in the following table. Some are defined 
as constants in the {\code W3IS2MD} module, others can be adjusted with the {\F SIS2} namelist. 

{\centering
\begin{tabular}{l  c c  c}
Parameters                   & Symbol        &  namelist parameter &  default values \\
\hline
Minimum floe size            & $D_{\min}$    & N. A. & 20~m \\
Initial floe size            & $D_{init}$    & N. A. & 1000~m \\
Ice fragility                & $\xi $        & N. A. & $0.9$ \\
Ice density                  & $\rho_{ice}$  & N. A. & 922.5~kg~m$^{-3}$ \\
Effective Young Modulus      & $Y^{*}$       & N. A. & 5.49~GPa \\
Poisson Coefficient          & $\nu$         & N. A. & 0.3 \\
Flexural strength            &  $\sigma_c$   & N. A. & 0.27~MPa \\
Flow law parameter           &$n$            & {\code IS2CREEPN} & 3\\
Flow law parameter           &$B$            & {\code IS2CREEPB} &  10$^7$~s$^{1/3}$\\
Elastic energy correction    & $C_R$         & {\code IS2BREAKE} & 1.0 \\
Relax. of disloc. compliance & $\delta D^{d}$& {\code IS2ANDISD} & $\Delta_d \Omega b^2/K$ (Pa$^{-1}$)\\
Dislocation density          & $\Delta_d$    & N. A. & 1.8$\times10^{9}$ (m$^{-2}$)\\
Restoring stress term        & $K$           & N. A. & 0.07 (Pa) \\
Orientation factor           & $\Omega$      & N. A. & $\pi^{-1}$ \\
Burgers vector               & $b$           & N. A. & 4.52$\times10^{-10}$ (m) \\
Drag term                    & $B$           & N. A. & $B_0 \exp(Q_v/(k_bT_K)$ (Pa.s)\\
 -                           & $B_0$         & N. A. & 1.205$\times10^{-9}$ (Pa.s)\\
 Boltzmann constant          & $k_B$         & N. A. & 8.617$\times10^{-5}$ (eV.K$^{-1}$) \\
Activation energy            & $Q_v$         & {\code IS2ANDISE} & 0.55 (eV) \\
 Peak broadening term        & $\alpha_{d}$  & N. A. & 0.54 \\
 Reduced variable            & $s$           & N. A. & $\log(\tau\omega)$ \\
 Relax. time of disloc.      & $\tau$        & N. A. & $B/K$ (s) \\
 Temperature of sea ice      & $T_K$         & N. A. & $268.15$ (K)\\
\hline
\end{tabular}
}



